!     ------------------
!     PROGRAMME ECLIPSOL  JMG 9904
!     ------------------
!
!     Jean Chapront,  Michelle Chapront-Touze, Gerard Francou.
!     Observatoire de Paris / DANOF / 9904
!
! --- Objet ------------------------------------------------------------
!
!     Recherche des eclipses centrales de Soleil entre -3000 et +3000
!     et determination de leurs lignes de centralite.
!
! --- Principales specifications ---------------------------------------
!
!     Les tables de la Lune utilisees sont issues de ELP2000-85, a semi-
!     analytical lunar ephemeris adequate for historical times,
!     Chapront Touze M., Chapront J., 1988, Astron. Astrophys., 190,342.
!
!     Les tables du Soleil utilisees sont celles de P. Bretagnon et de
!     J.L. Simon issues de la solution VSOP82 (Astron. Astrophys., 1982,
!     144,278) et de termes seculaires de la theorie J. Laskar (Astron.
!     Astrophys., 1986, 157, 59).
!
!     Le calcul des dates et lieux des eclipses est effectue a l'aide de
!     de representations polynomiales des coordonnees des corps (Lune et
!     Soleil) et des elements de Bessel.
!
!     La localisation des points de la ligne de centralite d'une eclipse
!     se calcule a partir des coordonnees apparentes de la Lune et du
!     Soleil rapportees dans le repere equatorial vrai de la date.
!
!     Le calcul de la duree du passage de l'ombre et de sa largeur ainsi
!     que la conversion entre dates calendaires et dates juliennes sont
!     effectues en utilisant des algorithmes de J. Meeus.
!
!     La difference entre Temps Terrestre TT et Temps Universel UT est
!     estimee de manieres differentes selon l'epoque concernee.
!     Avant 1955, a partir des tables de Stephenson et Morrison (1984).
!     Entre 1995 et 2000, a partir des tables TAI-UTC de l'IERS.
!     Apres 2000, a partir d'une extrapolation des valeurs 1800-2000.
!
! --- Fichiers ---------------------------------------------------------
!
!     tablun.dat   : Unite 11 (Lecture)   Tables de la Lune.
!     eclipsol.txt : Unite 12 (Ecriture)  Sauvegarde des resultats.
!
! --- Subroutines ------------------------------------------------------
!
!     LIRLUN : Lecture des tables de la Lune.
!     CALCUL : Calcul des eclipses centrales d'une annee.
!     RESULT : Affichage et sauvegarde des resultats.
!
!     SUBLUN : Substitution du temps dans les tables de la Lune.
!     SUBSOL : Substitution du temps dans les tables du Soleil.
!     POLYNO : Calcul des polynomes d'approximation des variables.
!     SUBPOL : Substitution du temps dans les polynomes.
!
!     LIGNE  : Coordonnees de la ligne de centralite de l'eclipse.
!     OMBRE  : Duree de l'ombre et largeur de la bande d'ombre.
!
!     EQUAPP : Calcul des coordonnees apparentes (Lune et Soleil).
!     BESSEL : Calcul des elements de Bessel de l'eclipse.
!     NUTAT  : Calcul de la nutation en longitude et en obliquite.
!     TSIDV  : Calcul du Temps Sideral vrai aa meridien origine.
!     DTTUT  : Difference Temps Terrestre (TT) - Temps Universel (UT).
!
!     DATEDJ : Conversion  Date calendaire ->  Date julienne.
!     DJDATE : Conversion  Date julienne   ->  Date calendaire.
!     CONTV  : Continuite d'une suite de valeurs en radians.
!     EFFECR : Effacement de l'ecran.
!
! --- Quantites transmises par common ----------------------------------
!
!     Tables de la Lune lues dans LIRLUN et utilisees dans SUBLUN.
!
!     common/lun/
!    .w(5),pa(7),
!    .pbc(250,3),cpbc(250,5,3),
!    .per(250,3),cper(250,2,3),
!    .pert(250,3),cpert(250,2,3),
!    .pert2(50,3),cpert2(50,2,3),
!    .aber(2,3),caber(2,2,3),
!    .kc(3),kp(3),kt(3),k2(3)
!
!     Tables des coefficients des polynomes Tchebychev calculees dans
!     POLYNO et utilisees dans SUBPOL.
!
!     common/pol/
!    .cal(10,3),cas(10,3),celb(10,9)
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      integer an
      character rep
      character*20 z,zb
!
      dimension resul(3,7,120),npt(3)
!
      parameter (d2000=2451545.d0,siecle=36525.d0)
!
! --- Initialisation ---------------------------------------------------
!
      nul1=11
      nul2=12
!
      open (nul1,file='tablun.dat',status='old')
      open (nul2,file='eclipsol.txt')
!
      write (zb,'(20x)')
!
! --- Lecture des tables de la Lune ------------------------------------
!
      write (*,"(//2x,'Lecture des tables de la Lune ... '//)")
!
      call LIRLUN (nul1)
!
! --- Debut du traitement ----------------------------------------------
!
      do
!
! --- Choisir l'annee --------------------------------------------------
!
         call EFFECR
         write (*,1000)
!
         write (*,1001)
         read (*,'(a)',iostat=nerr) z
         if (z(1:1).eq.' ') cycle
         read (z,'(i8)',iostat=nerr) an
!
         if (nerr.ne.0.or.an.lt.-3000.or.an.gt.3000) then
            write (*,1002)
            read (*,*)
            cycle
         endif
!
! --- Choisir l'echelle de temps (TT ou UT) ----------------------------
!
         do
            call EFFECR
            write (*,1000)
            write (*,1003) an
            read (*,'(i8)',iostat=nerr) it
            if (nerr.ne.0.or.it.lt.1.or.it.gt.2) cycle
            exit
         enddo
!
! --- Confirmation des choix -------------------------------------------
!
         write (*,1005)
         read (*,'(a)') rep
!
         if (rep.eq.'n'.or.rep.eq.'N') cycle
!
         write (*,1006)
!
! --- Eclipses centrales pour l'annee et l'echelle de temps choisies ---
!
         do
!
            ka=an
            call DATEDJ (1,1,ka,0,0,0,dj)
            tdeb=(dj-d2000)/siecle
!
            ka=an+1
            call DATEDJ (1,1,ka,0,0,0,dj)
            tfin=(dj-d2000)/siecle
!
            select case (it)
            case (1)
               dift=0.d0
            case (2)
               do
                  call EFFECR
                  write (*,1000)
                  call DATEDJ (1,7,an,0,0,0,t)
                  call DTTUT (t,dift)
                  write (*,1004) an,dift
                  write (z,'(20x)')
                  read (*,'(a)') z
                  if (z(1:1).eq.' ') then
                     exit
                  else
                     if (index(z,'.').eq.0) cycle
                     read (z,'(f20.5)',iostat=nerr) d
                     if (nerr.ne.0.or.d.lt.-5.d6.or.d.gt.5.d6) cycle
                     dift=d
                     exit
                  endif
               enddo
            endselect
!
            call CALCUL (tdeb,tfin,dift,necl,npt,resul)
            call RESULT (an,it,dift,nul2,necl,npt,resul)
!
            do
               call EFFECR
               write (*,1000)
               write (*,1007) an
               read (*,'(a)') rep
               if (rep.eq.' ') cycle
               if (rep.eq.'.') stop
               exit
            enddo

            if (rep.eq.'+'.or.rep.eq.'-') then
               if (rep.eq.'+') an=min0(an+1,+3000)
               if (rep.eq.'-') an=max0(an-1,-3000)
               cycle
            endif
!
            exit
!
         enddo
!
         write (*,1008)
         read (*,'(a)') rep
         if (rep.eq.'n'.or.rep.eq.'N') stop
!
! --- Fin du traitement ------------------------------------------------
!
      enddo
!
! --- Formats ----------------------------------------------------------
!
1000  format (2x,69('-')/
     .        2x,'ECLIPSES CENTRALES DE SOLEIL ENTRE -3000 ET +3000',
     .        2x,'(ECLIPSOL JMG9904)'/
     .        2x,69('-')/)
1001  format (/2x,'Annee (de -3000 a +3000) ?  ',$)
1002  format (/2x,'Erreur de reponse : Appuyez sur Entree ')
1003  format (/2x,'Annee : ',i5/
     .        /2x,'Temps Terrestre (1) ou Temps Universel (2) ?  ',$)
1004  format (/2x,'Annee : ',i5/
     .        /2x,'Valeur estimee de TT-UT  en secondes : ',f9.1/
     .        /2x,'Pour confirmer cette valeur, taper sur Entree'/
     .        /2x,'Sinon, Entrez la valeur de votre choix'
     .        /9x,'sans oublier le point decimal :    ',$)
1005  format (//2x,'Confirmez-vous vos choix (o/n) ? ',$)
1006  format (//2x,'Calcul en cours ... ')
1007  format (//2x,'Fin du calcul Annee ',i5/
     .        /2x,'Arret (.)',4x,'Suivante (+)',4x,'Precedente (-)',
     .         4x,'Autre Annee (*)',4x,'?  ',$)
1008  format (//2x,'Voulez-vous faire un autre calcul (o/n) ? ',$)
!
      end
!
!
!
      subroutine LIRLUN (nul)
! ----------------------------------------------------------------------
!
!     Ref : JMG 9904
!
!     Lecture des tables de la Lune.
!
!     Source : ELP2000-85, a semi-analytical lunar ephemeris for
!     for historical times, Chapront Touze M., Chapront J., 1988,
!     Astron. Astrophys., 190, 342.
!
! --- Entree ----------------------------------------------------------
!
!     nul  : Numero d'unite logique du fichier tablun.dat (entier).
!
! --- Common en sortie : lun -------------------------------------------
!
!     Longitude moyenne et precession.
!     Series du probleme central.
!     Series des perturbations (termes periodiques).
!     Series des perturbations (termes mixtes ordre 1).
!     Series des perturbations (termes mixtes ordre 2).
!     Termes de l'aberration.
!     Nombres de termes des series par variable.
!
! --- Remarque ---------------------------------------------------------
!
!     Cette subroutine lit les series du fichier et constitue les
!     tables de termes qui permettront la substitution du temps pour
!     calculer des coordonnees moyennes de la date :
!     Longitude et Latitude geocentriques et Distance Terre-Lune.
!     Elle fournit egalement les coefficients pour la correction de
!     l'aberration.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension ilu(4),ipla(13),phi(5)
      dimension ad0(4),ad1(4),ad2(4),ad3(4),ad4(4)
      dimension pa0(7),apla0(13),apla1(13)
      dimension unit(3)
!
      common/lun/
     .w(5),pa(7),
     .pbc(250,3),cpbc(250,5,3),
     .per(250,3),cper(250,2,3),
     .pert(250,3),cpert(250,2,3),
     .pert2(50,3),cpert2(50,2,3),
     .aber(2,3),caber(2,2,3),
     .kc(3),kp(3),kt(3),k2(3)
!
! --- Constantes -------------------------------------------------------
!
      parameter (ua=1.d0/1.49597870d08)
      parameter (pi=3.141592653589793d0,sdrad=pi/180.d0/3600.d0)
      parameter (d90=90.d0*3600.d0)
!
! --- Precession (") ---------------------------------------------------
!
      data pa0/
     .50287.275d0,110.5394d0,0.076d0,-0.2353d0,-0.00181d0,
     .0.000175d0,0.000013d0/
!
! --- Moyen mouvement des planetes (") ---------------------------------
!
      data apla1/
     . 538101628.68898d0, 210664136.43355d0,
     . 129597742.27580d0,  68905077.59284d0, 10925660.42861d0,
     .   4399609.65932d0,   1542481.19393d0,   786550.32074d0, 5*0.d0/
!
! --- Longitude moyenne de la Lune (") ---------------------------------
!
      w(1)=     785939.83019d0
      w(2)= 1732559343.34732d0
      w(3)=         -6.8098d0
      w(4)=         +0.006604d0
      w(5)=         -0.00003169d0
!
! --- Termes de la precession en "/siecle -----------------------------
!
      diz=1.d0
      do i=1,7
         diz=10*diz
         pa(i)=pa0(i)/diz
      enddo
!
! --- Arguments de Delaunay (") ----------------------------------------
!
      ad0(1)= 297.85019204d0*3600.d0
      ad0(2)= 357.52910254d0*3600.d0
      ad0(3)= 134.96339641d0*3600.d0
      ad0(4)=  93.27209240d0*3600.d0
!
      ad1(1)= 1602961601.0485d0
      ad1(2)=  129596581.0705d0
      ad1(3)= 1717915923.0161d0
      ad1(4)= 1739527262.9526d0
!
      ad2(1)=  -6.7896d0
      ad2(2)=  -0.5529d0
      ad2(3)= +31.4541d0
      ad2(4)= -13.1691d0
!
      ad3(1)= +0.006595d0
      ad3(2)= +0.000147d0
      ad3(3)= +0.051651d0
      ad3(4)= -0.001021d0
!
      ad4(1)= -0.00003184d0
      ad4(2)= +0.00000015d0
      ad4(3)= -0.00024470d0
      ad4(4)= +0.00000417d0
!
! --- Arguments des perturbations (") ----------------------------------
!
      apla0(1)= 252*3600.d0 +15*60.d0 + 3.25986d0
      apla0(2)= 181*3600.d0 +58*60.d0 +47.28305d0
      apla0(3)= 100*3600.d0 +27*60.d0 +59.22059d0
      apla0(4)= 355*3600.d0 +25*60.d0 +59.78866d0
      apla0(5)=  34*3600.d0 +21*60.d0 + 5.34212d0
      apla0(6)=  50*3600.d0 + 4*60.d0 +38.89694d0
      apla0(7)= 314*3600.d0 + 3*60.d0 +18.01841d0
      apla0(8)= 304*3600.d0 +20*60.d0 +55.19575d0
      apla0(9)= w(1)
      apla1(9)= w(2)+pa(1)
      do i=10,13
         apla0(i)=ad0(i-9)
         apla1(i)=ad1(i-9)
      enddo
!
! --- Constantes pour la conversion d'unite des 3 coordonnees ---------
!
      unit(1)=sdrad
      unit(2)=sdrad
      unit(3)=ua
!
! --- Lecture des series du probleme central---------------------------
!
      do i=1,3
         coef=unit(i)
         read (nul,1001) kc(i)
         do j=1,kc(i)
            read (nul,1002) (ilu(k),k=1,4),dp
            pbc(j,i)=dp*coef
            do k=1,5
               phi(k)=0.d0
            enddo
            do k=1,4
               phi(1)=phi(1)+ilu(k)*ad0(k)
               phi(2)=phi(2)+ilu(k)*ad1(k)
               phi(3)=phi(3)+ilu(k)*ad2(k)
               phi(4)=phi(4)+ilu(k)*ad3(k)
               phi(5)=phi(5)+ilu(k)*ad4(k)
            enddo
            if (i.eq.3) then
               phi(1)=phi(1)+d90
            endif
            do k=1,5
               cpbc(j,k,i)=phi(k)*sdrad
            enddo
         enddo
      enddo
!
! --- Lecture des series des perturbations (termes periodiques) -------
!
      do i=1,3
         coef=unit(i)
         read (nul,1001) kp(i)
         do j=1,kp(i)
            read (nul,1003) (ipla(k),k=1,13),pha,dp
            per(j,i)=dp*coef
            phi0=0.d0
            phi1=0.d0
            do k=1,13
               phi0=phi0+apla0(k)*ipla(k)
               phi1=phi1+apla1(k)*ipla(k)
            enddo
            phi0=phi0+pha*3600.d0
            cper(j,1,i)=phi0*sdrad
            cper(j,2,i)=phi1*sdrad
         enddo
      enddo
!
! --- Lecture des series des perturbations (termes mixtes ordre 1) ----
!
      do i=1,3
         coef=unit(i)
         read (nul,1001) kt(i)
         do j=1,kt(i)
            read (nul,1003) (ipla(k),k=1,13),pha,dp
            pert(j,i)=dp*coef
            phi0=0.d0
            phi1=0.d0
            do k=1,13
               phi0=phi0+apla0(k)*ipla(k)
               phi1=phi1+apla1(k)*ipla(k)
            enddo
            phi0=phi0+pha*3600.d0
            cpert(j,1,i)=phi0*sdrad
            cpert(j,2,i)=phi1*sdrad
         enddo
      enddo
!
! --- Lecture des series des perturbations (termes mixtes ordre 2) ----
!
      do i=1,3
         coef=unit(i)
         read (nul,1001) k2(i)
         do j=1,k2(i)
            read (nul,1003) (ipla(k),k=1,13),pha,dp
            pert2(j,i)=dp*coef
            phi0=0.d0
            phi1=0.d0
            do k=1,13
               phi0=phi0+apla0(k)*ipla(k)
               phi1=phi1+apla1(k)*ipla(k)
            enddo
            phi0=phi0+pha*3600.d0
            cpert2(j,1,i)=phi0*sdrad
            cpert2(j,2,i)=phi1*sdrad
         enddo
      enddo
!
! --- Lecture des termes de l'aberration ------------------------------
!
      do i=1,3
      do j=1,2
         aber(j,i)=0.d0
         do k=1,2
            caber(j,k,i)=0.d0
         enddo
      enddo
      enddo
!
      aber(1,1)=-0.703d0*sdrad
      caber(1,1,1)=d90*sdrad
!
      aber(2,1)=-0.038d0*sdrad
      caber(2,1,1)=(ad0(3)+d90)*sdrad
      caber(2,2,1) =ad1(3)*sdrad
!
      aber(1,2)=-0.063d0*sdrad
      caber(1,1,2)=(ad0(4)+d90)*sdrad
      caber(1,2,2)=ad1(4)*sdrad
!
      aber(1,3)=-0.068*ua
      caber(1,1,3)=ad0(3)*sdrad
      caber(1,2,3)=ad1(3)*sdrad
!
      aber(2,3)=-0.23*ua
      caber(2,1,3)=2*ad0(1)*sdrad
      caber(2,2,3)=2*ad1(1)*sdrad
!
      return
!
! --- Formats ---------------------------------------------------------
!
1001  format (i4)
1002  format (3x,4i3,2x,f13.5)
1003  format (3x,13i3,2x,f9.5,2x,f9.5)
!
      end
!
!
!
      subroutine CALCUL (tdeb,tfin,dift,necl,npt,resul)
! ----------------------------------------------------------------------
!
!     Ref : JMG 9904
!
!     Calcul des eclipses centrales de Soleil entre deux dates.
!     Parametres de la ligne de centralite de chaque eclipse.
!
! --- Entree -----------------------------------------------------------
!
!     tdeb : Date du debut de l'annee en siecle julien depuis J2000
!            (reel double).
!     tfin : Date de la fin de l'annee en siecle julien depuis J2000
!            (reel double).
!
!     dift : Difference TT-UT en secondes (reel double).
!
! --- Sortie -----------------------------------------------------------
!
!     necl  : Nombre d'eclipses (entier, max:3).
!     npt   : Nombre de points de la ligne (entier,max:120).
!     resul : Table des resultats par eclipse (reel double).
!
!     Pour chaque eclipse centrale, la table resul contient les donnees
!     concernant les points de la ligne de centralite :
!
!     1 : Date en siecle julien depuis J2000 (reel dp).
!     2 : Longitude du lieu en degres (reel double).
!     3 : Latitude du lieu en degres (reel double).
!     4 : Rayon du cone d'ombre en km (reel double).
!     5 : Largeur de la bande de centralite en km (reel double).
!     6 : Duree de l'eclipse en minutes (reel double).
!     7 : Hauteur du Soleil en degres (reel double).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      double precision lonl,latl,lons,lats,ls0,ls1,latmin
!
      dimension resul(3,7,120),npt(3)
      dimension rm(7,100),rd(7,10),rf(7,10)
      dimension v1(9),v2(9)
!
! --- Constantes -------------------------------------------------------
!
      parameter (d2000=2451545.d0,siecle=36525.d0)
      parameter (pi=3.141592653589793d0,dpi=2.d0*pi)
      parameter (dgrad=pi/180.d0,sdrad=dgrad/3600.d0)
!
! --- Precision pour la determination de la Nouvelle Lune -------------
!
!     Precision 1 : 1 heure  en siecle julien.
!     Precision 2 : 1 minute en siecle julien.
      parameter (prec1=1.d0/24.d0/siecle,prec2=prec1/60.d0)
!
! --- Valeur limite de la latitude de la Lune --------------------------
!
!     Limite : 1.038 degre en radian.
      parameter (latmin=1.038d0*dgrad)
!
! --- Ligne de centralite date par date --------------------------------
!
!     Nombre de points : 96.
!     Pas 1 : 5 minutes en siecle julien.
!     Pas 2 : 1 minute  en siecle julien.
      parameter (npas=96,pas=5.d0)
      parameter (pas1=pas/60.d0/24.d0/36525.d0,pas2=pas1/pas)
!
! --- Parametres des polynomes de Tchebychev ---------------------------
!
!     Nombre de coefficients (ncf) : 4.
!     Intervalle de validite (delta) : 8 heures en siecle julien.
      parameter (ncf=4,delta=8.d0/24.d0/36525.d0)
!
! --- Initialisation ---------------------------------------------------
!
      necl=0
      dte=dift/86400.d0/siecle
!
! --- Longitudes moyennes de la Lune et du Soleil ----------------------
!
      al0=785939.83019d0*sdrad
      al1=(1732559343.34732d0+5028.7275d0)*sdrad
      as0=4.9353929d0
      as1=62833.1961680d-2-0.3691d0*sdrad
      ls0=al0-as0
      ls1=al1-as1
      akl=(ls0+ls1*tdeb)/dpi
      kl=int(akl)-1
!
! --- Debut du traitement ----------------------------------------------
!
      do
!
! --- Recherche des Nouvelles Lunes (precision 1) ----------------------
!
         t=(kl*dpi-ls0)/ls1
         dt=1.d0
!
         do while (abs(dt).gt.prec1)
            call SUBLUN (t,lonl,latl,r,1)
            call SUBSOL (t,lons,lats,r)
            dt=(kl*dpi-lonl+lons)/ls1
            t=t+dt
         enddo
!
! --- Nouvelle Lune hors limite ----------------------------------------
!
         if (t.gt.tfin) return
!
         if (t.lt.tdeb) then
            kl=kl+1
            cycle
         endif
!
! --- Recherche des Nouvelles Lunes (precision 2) ----------------------
!
         dt=1.d0
!
         do while (abs(dt).gt.prec2)
            call SUBLUN (t,lonl,latl,r,1)
            call SUBSOL (t,lons,lats,r)
            dt=(kl*dpi-lonl+lons)/ls1
            t=t+dt
         enddo
!
! --- Test sur la latitude de la Lune ----------------------------------
!
         call SUBLUN (t,lonl,latl,r,2)
!
         if (dabs(latl).gt.latmin) then
            kl=kl+1
            cycle
         endif
!
! --- Calcul des polynomes representant les variables du calcul --------
!
         tu=t-delta/2.d0-dte
         t0=aint(tu/pas1)*pas1+dte
!
         call POLYNO (t0,delta,ncf)
!
         nm=0
         nd=0
         nf=0
!
! --- Ligne de centralite ----------------------------------------------
!
         kle=1
!
         do it=1,npas
            te=t0+(it-1)*pas1
            call SUBPOL (1,te,t0,delta,ncf,v1,v2)
            call LIGNE (te,v1,v2,glon,glat,rayc,haut,ierr)
            if (ierr.eq.0) then
               nm=nm+1
               call SUBPOL (2,te,t0,delta,ncf,v1,v2)
               call OMBRE (v1,v2,rayc,duree,bande)
               rm(1,nm)=te
               rm(2,nm)=glon
               rm(3,nm)=glat
               rm(4,nm)=rayc
               rm(5,nm)=duree
               rm(6,nm)=bande
               rm(7,nm)=haut
               if (kle.eq.1) then
                  ted=te
                  kle=2
               endif
            else
               if (kle.eq.2) then
                  tef=te-pas1
                  kle=3
                  exit
               endif
            endif
         enddo
!
         if (kle.eq.1) then
            kl=kl+1
            cycle
         endif
!
         if (kle.eq.2) tef=te
!
! --- Debut de la ligne de centralite ----------------------------------
!
         do it=1,9
            te=ted-it*pas2
            if (te.lt.t0) exit
            call SUBPOL (1,te,t0,delta,ncf,v1,v2)
            call LIGNE (te,v1,v2,glon,glat,rayc,haut,ierr)
            if (ierr.ne.0) exit
            call SUBPOL (2,te,t0,delta,ncf,v1,v2)
            call OMBRE (v1,v2,rayc,duree,bande)
            nd=nd+1
            rd(1,nd)=te
            rd(2,nd)=glon
            rd(3,nd)=glat
            rd(4,nd)=rayc
            rd(5,nd)=duree
            rd(6,nd)=bande
            rd(7,nd)=haut
         enddo
!
! --- Fin de la ligne de centralite ------------------------------------
!
         do it=1,9
            te=tef+it*pas2
            if (te.gt.t0+delta) exit
            call SUBPOL (1,te,t0,delta,ncf,v1,v2)
            call LIGNE (te,v1,v2,glon,glat,rayc,haut,ierr)
            if (ierr.ne.0) exit
            call SUBPOL (2,te,t0,delta,ncf,v1,v2)
            call OMBRE (v1,v2,rayc,duree,bande)
            nf=nf+1
            rf(1,nf)=te
            rf(2,nf)=glon
            rf(3,nf)=glat
            rf(4,nf)=rayc
            rf(5,nf)=duree
            rf(6,nf)=bande
            rf(7,nf)=haut
         enddo
!
! --- Rangement des resultats ------------------------------------------
!
         necl=necl+1
         npt(necl)=nm+nd+nf
!
         do i=1,nd
            j=nd-i+1
            do k=1,7
               resul(necl,k,i)=rd(k,j)
            enddo
         enddo
!
         do i=1,nm
            j=i+nd
            do k=1,7
               resul(necl,k,j)=rm(k,i)
            enddo
         enddo
!
         do i=1,nf
            j=i+nd+nm
            do k=1,7
               resul(necl,k,j)=rf(k,i)
            enddo
         enddo
!
         kl=kl+1
!
! --- Fin du traitement ------------------------------------------------
!
      enddo
!
      return
      end
!
!
!
      subroutine RESULT (an,it,dift,nul,necl,npt,resul)
! ----------------------------------------------------------------------
!
!     Ref : JMG 9904
!
!     Affichage et Sauvegarde des resultats.
!
! --- Entree -----------------------------------------------------------
!
!     an    : Annee (entier).
!     it    : Indice de l'echelle de temps (entier, 1:TT 2:UT).
!     dift  : Difference TT-UT en secondes (reel double).
!     nul   : Unite Logique du fichier eclipsol.txt (entier).
!
!     necl  : Nombre d'eclipses (entier, max:3).
!     npt   : Nombre de points de la ligne  (entier, max:120).
!     resul : Table des resultats par eclipse (reel double).
!
!     Pour chaque eclipse centrale, la table resul contient les donnees
!     concernant les points de la ligne de centralite :
!
!     1 : Date en siecle julien depuis J2000 (reel dp).
!     2 : Longitude du lieu en degres (reel double).
!     3 : Latitude du lieu en degres (reel double).
!     4 : Rayon du cone d'ombre en km (reel double).
!     5 : Duree de l'eclipse en minutes (reel double).
!     6 : Largeur de la bande de centralite en km (reel double).
!     7 : Hauteur du Soleil en degres (reel double).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character alo,ala
      character*2 ax(2)
      character*4 adur,alar,ahau
      character*9 mois(12),type(-1:1)
      character*16 adift
      character*25 at(2),am(2)
      integer an
      logical deb
!
      dimension resul(3,7,120),npt(3)
!
      parameter (d2000=2451545.d0,siecle=36525.d0)
!
! --- Initialisation ---------------------------------------------------
!
      data mois/
     .'JANVIER','FEVRIER','MARS','AVRIL','MAI','JUIN',
     .'JUILLET','AOUT','SEPTEMBRE','OCTOBRE','NOVEMBRE','DECEMBRE'/
!
      data type/'ANNULAIRE','MIXTE','TOTALE'/
!
      data at/'Temps Terrestre (TT)','Temps Universel (UT)'/
!
      data am/'Meridien des Ephemerides','Meridien de Greenwich'/
!
      data ax/'TT','UT'/
!
      data deb/.true./
!
      dte=dift/86400.d0
!
      if (deb) then
         write (nul,1000)
         deb=.false.
      endif
!
! --- Absence d'eclipse centrale ---------------------------------------
!
      if (necl.eq.0) then
         call EFFECR
         write (*,1001)
         write (*,1002) an
         write (nul,1002) an
         write (*,1005)
         read (*,*)
         return
      endif
!
! --- Boucle Eclipse : Debut -------------------------------------------
!
      do i=1,necl
!
! --- Date de l'eclipse ------------------------------------------------
!
         n=npt(i)/2
         te=resul(i,1,n)*siecle+d2000
         tu=te-dte
         dj=aint(tu*1440.d0+0.5)/1440.d0
         call DJDATE (dj,0,ij0,im0,ia0,nh,nm,ns)
!
! --- Type de l'eclipse ------------------------------------------------
!
         ityp=1
         if (resul(i,4,1).lt.0.d0) ityp=-1
         do j=1,npt(i)
            k=ityp*resul(i,4,j)
            if (k.lt.0) ityp=0
         enddo
!
! --- Difference TT-UT -------------------------------------------------
!
         if (it.eq.1) then
            adift=' '
         else
            write (adift,"('TT-UT = ',f7.1,'s')") dift
            if (adift(12:14).eq.'  .') adift(12:14)=' 0.'
            if (adift(12:14).eq.' -.') adift(12:14)='-0.'
         endif
!
         write (nul,'(/)')
         write (nul,1003) an,ij0,mois(im0),type(ityp),adift,
     .                    at(it),am(it),ax(it)
!
! --- Ligne de centralite point par point ------------------------------
!
         do j=1,npt(i)
!
            te=resul(i,1,j)*siecle+d2000
            glon=resul(i,2,j)
            glat=resul(i,3,j)
            rayc=resul(i,4,j)
            duree=resul(i,5,j)
            bande=resul(i,6,j)
            haut=resul(i,7,j)
!
            if (mod(j,10).eq.1) then
               call EFFECR
               write (*,1001)
               write (*,1003) an,ij0,mois(im0),type(ityp),adift,
     .                        at(it),am(it),ax(it)
            endif
!
! --- Date du point de centralite --------------------------------------
!
            tu=te-dte
            dj=aint(tu*1440.d0+0.5)/1440.d0
            call DJDATE (dj,2,ij,im,ia,nh,nm,ns)
!
! --- Longitude du lieu ------------------------------------------------
!
            x=glon-dte*360.d0*1.0027379d0
            x=mod(x,360.d0)
            if (x.gt.+180.d0) x=x-360.d0
            if (x.le.-180.d0) x=x+360.d0
            glon=x
!
            k=abs(glon)*60.d0+0.5d0
            lod=k/60
            lom=k-60*lod
            alo=' '
            if (glon.gt.0) alo='O'
            if (glon.lt.0) alo='E'
!
! --- Latitude du lieu -------------------------------------------------
!
            k=abs(glat)*60.d0+0.5d0
            lad=k/60
            lam=k-60*lad
            ala=' '
            if (glat.gt.0) ala='N'
            if (glat.lt.0) ala='S'
!
! --- Rayon de la section du cone d'ombre ------------------------------
!
            kr=abs(rayc)+0.5d0
!
! --- Duree de l'eclipse -----------------------------------------------
!
            write (adur,'(f4.1)') duree
            if (adur(2:2).eq.' ') adur(2:2)='0'
!
! --- Largeur de la bande de centralite --------------------------------
!
            if (haut.ge.30.05d0) then
               b=bande+2.5d0
               l=int(b/5)*5
               write (alar,'(i4)') l
            else
               alar='  - '
            endif
!
! --- Hauteur du Soleil -----------------------------------------------
!
            write (ahau,'(f4.1)') haut
            if (ahau(2:2).eq.' ') ahau(2:2)='0'
!
! --- Affichage et Sauvegarde des quantites liees au point -------------
!
            write (nul,1004) ij,im,ia,nh,nm,
     .                       lad,lam,ala,lod,lom,alo,kr,adur,alar,ahau
            write (*,1004)   ij,im,ia,nh,nm,
     .                       lad,lam,ala,lod,lom,alo,kr,adur,alar,ahau
!
            if (mod(j,10).eq.0) then
               write (*,1005)
               read (*,*)
            endif
!
         enddo
!
         if (mod(npt(i),10).ne.0) then
            write (*,1005)
            read (*,*)
         endif
!
! --- Boucle Eclipse : Fin ---------------------------------------------
!
      enddo
!
      return
!
! --- Formats ----------------------------------------------------------
!
1000  format (2x,69('-')/
     .        2x,'ECLIPSES CENTRALES DE SOLEIL ENTRE -3000 ET +3000',
     .        2x,'(ECLIPSOL JMG9904)'/
     .        2x,69('-')//
     .        2x,'Dates et heures des points de la ligne de',
     .        1x,'centralite c''est-a-dire des'/
     .        2x,'points de contact de l''axe du cone d''ombre avec',
     .        1x,'la surface terrestre.'//
     .        2x,'Lat : Latitude geographique du point par rapport',
     .        1x,'au meridien origine'/
     .        2x,'Lon : Longitude geographique du point par rapport',
     .        1x,'au meridien origine'/
     .        2x,'Ray : Rayon de la section normale du cone d''ombre'/
     .        2x,'Dur : Duree de la phase centrale de l''eclipse'/
     .        2x,'Lar : Largeur estimee de la bande d''ombre'/
     .        2x,'Hau : Hauteur du Soleil'//
     .        2x,'Echelle de temps : Temps Terrestre (TT)      / ',
     .        1x,'Temps Universel (UT)'/
     .        2x,'Meridien origine : Meridien des Ephemerides  / ',
     .        1x,'Meridien de Greenwich'/)
1001  format (2x,69('-')/
     .        2x,'ECLIPSES CENTRALES DE SOLEIL ENTRE -3000 ET +3000',
     .        2x,'(ECLIPSOL JMG9904)'/
     .        2x,69('-'))
1002  format (/2x,'ANNEE ',i5,4x,'Absence d''eclipse centrale')
1003  format (/2x,'ANNEE ',i5,4x,i2.2,1x,a9,3x,'ECLIPSE ',a9,6x,a16/
     .        /2x,'Echelle de temps : ',a25
     .        /2x,'Meridien origine : ',a25/
     .        /2x,'Date et Heure (',a2,')',8x,'Lat',8x,'Lon',7x,'Ray',
     .         6x,'Dur',6x,'Lar',5x,'Hau'/)
1004  format (2x,i2.2,'/',i2.2,'/',i5,2x,i2.2,'h',i2.2,'m',
     .        5x,i2,'d',i2.2,'''',a1,
     .        3x,i3,'d',i2.2,'''',a1,
     .        4x,i3,' km',2x,a4,' m',2x,a4,' km',3x,a4,'d')
1005  format (/2x,'Appuyez sur Entree')
!
      end
!
!
!
      subroutine SUBLUN (t,long,lat,r,k)
! ----------------------------------------------------------------------
!
!     Ref : JMG 9904
!
!     Substitution du temps dans les tables de la Lune.
!
!     Les coordonnees calculees sont la longitude, la latidude et la
!     distance a la Terre. Elles sont rapportees au repere ecliptique
!     moyen de la date et corrigees de l'aberration.
!
!     Source : tables de la Lune de M. Chapront-Touze et J. Chapront,
!     ELP2000-85, a semi-analytical lunar ephemeris for historical times
!     (Astron Astrophys, 1988,190,342).
!
! --- Entree -----------------------------------------------------------
!
!     t : Date TT en siecle julien depuis J2000 (reel double).
!
!     k : Indice des coordonnees (entier).
!         1 : Longitude.
!         2 : Longitude et Latitude.
!         3 : Longitude, Latitude et Distance Terre-Lune.
!
! --- Common en entree : lun -------------------------------------------
!
!     Longitude moyenne et precession.
!     Series du probleme central.
!     Series des perturbations (termes periodiques).
!     Series des perturbations (termes mixtes ordre 1).
!     Series des perturbations (termes mixtes ordre 2).
!     Termes de l'aberration.
!     Nombres de termes des series par variable.
!
! --- Sortie -----------------------------------------------------------
!
!     long : Longitude de la Lune en radian (reel double).
!     lat  : Latitude de la Lune en radian (reel double).
!     r    : Distance Terre-Lune en ua (reel double).
!
! --- Remarque ---------------------------------------------------------
!
!     Les termes des series du probleme central et des perturbations
!     sont fournis par common (voir subroutine LIRLUN).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      double precision long,lat
!
      dimension c(3),ct(3),ct2(3)
!
      parameter (pi=3.141592653589793d0,sdrad=pi/180.d0/3600.d0)
!
      common/lun/
     .w(5),pa(7),
     .pbc(250,3),cpbc(250,5,3),
     .per(250,3),cper(250,2,3),
     .pert(250,3),cpert(250,2,3),
     .pert2(50,3),cpert2(50,2,3),
     .aber(2,3),caber(2,2,3),
     .kc(3),kp(3),kt(3),k2(3)
!
! --- Substitution du temps dans les series par variable ---------------
!
      do i=1,k
!
! ------ Series du probleme central -----------------------------------
!
         c(i)=0.d0
!
         do j=1,kc(i)
            phi=cpbc(j,1,i)+t*(cpbc(j,2,i)+t*(cpbc(j,3,i)
     .         +t*(cpbc(j,4,i)+t*cpbc(j,5,i))))
            s=sin(phi)
            c(i)=c(i)+s*pbc(j,i)
         enddo
!
! ------ Series des perturbations (termes periodiques) ----------------
!
         do j=1,kp(i)
            phi=cper(j,1,i)+t*cper(j,2,i)
            s=sin(phi)
            c(i)=c(i)+s*per(j,i)
         enddo
!
! ------ Series des perturbations (termes mixtes ordre 1) -------------
!
         ct(i)=0.d0
!
         do j=1,kt(i)
            phi=cpert(j,1,i)+t*cpert(j,2,i)
            s=sin(phi)
            ct(i)=ct(i)+s*pert(j,i)
         enddo
!
! ------ Series des perturbations (termes mixtes ordre 2) -------------
!
         ct2(i)=0.d0
!
         do j=1,k2(i)
            phi=cpert2(j,1,i)+t*cpert2(j,2,i)
            s=sin(phi)
            ct2(i)=ct2(i)+s*pert2(j,i)
         enddo
!
         c(i)=c(i)+t*(ct(i)+t*ct2(i))
!
! ------ Correction de l'aberration -----------------------------------
!
         do j=1,2
            phi=caber(j,1,i)+t*caber(j,2,i)
            s=sin(phi)
            c(i)=c(i)+s*aber(j,i)
         enddo
!
      enddo
!
! --- Longitude moyenne et precession ----------------------------------
!
      al=0.d0
      do i=5,2,-1
         al=t*(al+w(i))
      enddo
!
      ap=0.d0
      do i=7,1,-1
         ap=t*(ap+pa(i))
      enddo
!
      long=c(1)+(w(1)+al+ap)*sdrad
      lat=c(2)
      r=c(3)
!
      return
      end
!
!
!
      subroutine SUBSOL (t,long,lat,r)
! ----------------------------------------------------------------------
!
!     Ref : JMG 9904
!
!     Substitution du temps dans les tables du Soleil.
!
!     Les coordonnees calculees sont la longitude, la latidude et la
!     distance a la Terre. Elles sont rapportees au repere ecliptique
!     moyen de la date et corrigees de l'aberration.
!
!     Source : tables du Soleil de P. Bretagnon et J.L. Simon issues de
!     la solution VSOP82 (Astron Astrophys, 1982,144,278) et de termes
!     seculaires de la theorie J. Laskar (Astron Astrophys, 1986,157,59).
!
! --- Entree -----------------------------------------------------------
!
!     t : Date TT en siecle julien depuis J2000 (reel double).
!
! --- Sortie -----------------------------------------------------------
!
!     long : Longitude du Soleil en radian (reel double).
!     lat  : Latitude du Soleil en radian (reel double).
!     r    : Distance Terre-Soleil en ua (reel double).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      double precision long,lat
!
      dimension il(50),al(50),bl(50),ir(50)
!
      parameter (pi=3.141592653589793d0,sdrad=pi/180.d0/3600.d0)
!
! --- Initialisation ---------------------------------------------------
!
      data il/
     . 403406, 195207, 119433, 112392,   3891,   2819,   1721,      0,
     .    660,    350,    334,    314,    268,    242,    234,    158,
     .    132,    129,    114,     99,     93,     86,     78,     72,
     .     68,     64,     46,     38,     37,     32,     29,     28,
     .     27,     27,     25,     24,     21,     21,     20,     18,
     .     17,     14,     13,     13,     13,     12,     10,     10,
     .     10,     10/
!
      data al/
     .  4.721964d0,   5.937458d0,  1.115589d0, 5.781616d0, 5.5474d0,
     .  1.5120d0 ,4.1897d0 ,1.163d0  ,5.415d0  ,4.315d0  ,4.553d0  ,
     .  5.198d0  ,5.989d0  ,2.911d0  ,1.423d0  ,0.061d0  ,2.317d0  ,
     .  3.193d0  ,2.828d0  ,0.52d0   ,4.65d0   ,4.35d0   ,2.75d0   ,
     .  4.50d0   ,3.23d0   ,1.22d0   ,0.14d0   ,3.44d0   ,4.37d0   ,
     .  1.14d0   ,2.84d0   ,5.96d0   ,5.09d0   ,1.72d0   ,2.56d0   ,
     .  1.92d0   ,0.09d0   ,5.98d0   ,4.03d0   ,4.27d0   ,0.79d0   ,
     .  4.24d0   ,2.01d0   ,2.65d0   ,4.98d0   ,0.93d0   ,2.21d0   ,
     .  3.59d0   ,1.50d0   ,2.55d0   /
!
      data bl/
     .   1.621043d0, 62830.348067d0, 62830.821524d0, 62829.634302d0,
     .125660.5691d0,125660.9845d0  , 62832.4766d0  ,     0.813d0   ,
     . 125659.310d0, 57533.850d0,   -33.931d0,777137.715d0, 78604.191d0,
     .      5.412d0, 39302.098d0,   -34.861d0,115067.698d0, 15774.337d0,
     .   5296.670d0, 58849.27d0 ,  5296.11d0 , -3980.70d0 , 52237.69d0 ,
     .  55076.47d0 ,   261.08d0 , 15773.85d0 ,188491.03d0 , -7756.55d0 ,
     .    264.89d0 ,117906.27d0 , 55075.75d0 , -7961.39d0 ,188489.81d0 ,
     .   2132.19d0 ,109771.03d0 , 54868.56d0 , 25443.93d0 ,-55731.43d0 ,
     .  60697.74d0 ,  2132.79d0 ,109771.63d0 , -7752.82d0 ,188491.91d0 ,
     .    207.81d0 , 29424.63d0 ,    -7.99d0 , 46941.14d0 ,   -68.29d0 ,
     .  21463.25d0 ,157208.40d0 /
!
      data ir/
     .    0,   -97597,   -59715,   -56188,    -1556,    -1126,     -861,
     .  941,-264,-163,   0, 309,-158,   0, -54,   0, -93, -20,   0, -47,
     .    0,   0, -33, -32,   0, -10, -16,   0,   0, -24, -13,   0,  -9,
     .    0, -17, -11,   0,  31, -10,   0, -12,   0,  -5,   0,   0,   0,
     .    0,   0,   0,  -9/
!
! --- Substitution du temps --------------------------------------------
!
      u=t/100.d0
      cl=0.d0
      cr=0.d0
!
      do i=1,50
         a=al(i)+u*bl(i)
         s=sin(a)
         c=cos(a)
         cl=cl+s*il(i)
         cr=cr+c*ir(i)
      enddo
!
      correc=(-0.3691d0+(-0.006606d0-0.000001d0*t)*t)*t
!
      long=4.9353929d0+62833.1961680d0*u+cl*1.d-7+correc*sdrad
!
      a=3.10d0+62830.14d0*u
      c=cos(a)
      long=long+1.d-7*(-993+17*c)
      lat=0.d0
      r=1.0001026d0+1.d-7*cr
!
      return
      end
!
!
!
      subroutine POLYNO (t0,dt,nc)
! ----------------------------------------------------------------------
!
!     Ref : JMG 9904
!
!     Calcul des coefficients des polynomes de Tchebychev.
!
! --- Entree -----------------------------------------------------------
!
!     t0 : Date du debut de l'intervalle en siecle julien depuis J2000.
!          (reel double).
!     dt : Taille de l'intervalle en siecle julien (reel double).
!     nc : Nombre de coefficients (entier, 3<=nc<=10).
!
! --- Common en sortie : pol -------------------------------------------
!
!     Tables des coefficients de polynomes de Tchebychev (reel double).
!     Le nombre maximum de coefficients par variables est : 10.
!
!     1/ Coefficients des coordonnees rectangulaires equatoriales X,Y,Z
!        de la Lune et du Soleil (reel double) :
!        cal(10,1:3) : X,Y,Z de la Lune en ua.
!        cas(10,1:3) : X,Y,Z du Soleil en ua.
!
!     2/ Coefficients des elements de Bessel.
!        Coordonnees horaires de l'axe passant par le centre de la Terre
!        et parallele a l'axe du cone d'ombre.
!        celb(10,1) : Angle Horaire en radian.
!        celb(10,2) : Declinaison en radian.
!        Coefficients des coordonnees rectangulaires de la Lune dans le
!        systeme de coordonnees fondamental de l'eclipse.
!        celb(10,3) : X en rayon terrestre.
!        celb(10,4) : Y en rayon terrestre.
!        celb(10,5) : Z en rayon terrestre.
!        Coefficients des demi-angles au sommet des cones.
!        celb(10,6) : Demi-angle du cone de penombre en radian.
!        celb(10,7) : Demi-angle du cone d'ombre en radian.
!        Coefficients des rayons des Sections des cones avec le plan
!        fondamental de l'eclipse.
!        celb(10,8) : Rayon du cone de penombre en rayon terrestre.
!        celb(10,9) : Rayon du cone d'ombre en rayon terrestre.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      double precision lonl,lons,latl,lats
!
      dimension x(10),tch(10,10)
      dimension val(3,10),vas(3,10),velb(9,10),ah(10)
!
      parameter (pi=3.141592653589793d0,sdrad=pi/180.d0/3600.d0)
!
      common/pol/
     .cal(10,3),cas(10,3),celb(10,9)
!
! --- Initialisation ---------------------------------------------------
!
      do i=1,nc
         a=(2.d0*i-1.d0)*pi/(2.d0*nc)
         x(i)=cos(a)
         tch(1,i)=1.d0
         tch(2,i)=x(i)
         do j=3,nc
            tch(j,i)=2*x(i)*tch(j-1,i)-tch(j-2,i)
         enddo
      enddo
!
! --- Verification du nombre de coefficients ---------------------------
!
      if (nc.lt.3.or.nc.gt.10) stop ' Erreur Polynomes de Tchebychev'
!
! --- Calcul des valeurs des variables pour les racines des polynomes ---
!
      do i=1,nc
         t=t0+(x(i)+1.d0)*dt/2.d0
         call NUTAT (t,psi,es,ec)
         call TSIDV (t,tsv)
         call SUBLUN (t,lonl,latl,rl,3)
         call SUBSOL (t,lons,lats,rs)
         call EQUAPP (lonl,latl,rl,psi,es,ec,val(1,i))
         call EQUAPP (lons,lats,rs,psi,es,ec,vas(1,i))
         call BESSEL (val(1,i),vas(1,i),tsv,velb(1,i))
      enddo
!
! --- Calcul des coefficients (1) --------------------------------------
!
      do k=1,3
         do j=1,nc
            a1=0.d0
            a2=0.d0
            do i=1,nc
               a1=a1+val(k,i)*tch(j,i)
               a2=a2+vas(k,i)*tch(j,i)
            enddo
            cal(j,k)=a1/2.d0
            cas(j,k)=a2/2.d0
         enddo
         cal(1,k)=cal(1,k)/2.d0
         cas(1,k)=cas(1,k)/2.d0
      enddo
!
! --- Calcul des coefficients (2) --------------------------------------
!
      do i=1,nc
         ah(i)=velb(1,i)
      enddo
      call CONTV (ah,nc,0.d0)
      do i=1,nc
         velb(1,i)=ah(i)
      enddo
!
      do k=1,9
         do j=1,nc
            a1=0.d0
            do i=1,nc
               a1=a1+velb(k,i)*tch(j,i)
            enddo
            celb(j,k)=a1/2.d0
         enddo
         celb(1,k)=celb(1,k)/2.d0
      enddo
!
      return
      end
!
!
!
      subroutine SUBPOL (iv,t,t0,dt,nc,v1,v2)
! ----------------------------------------------------------------------
!
!     Ref : JMG 9904
!
!     Substitution du temps dans les polynomes de Tchebychev.
!
! --- Entree -----------------------------------------------------------
!
!     iv : Indice des variables (entier).
!          iv=1 : Coordonnees equatoriales de la Lune et du Soleil.
!          iv=2 : Elements de Bessel de l'eclipse.
!
!     t  : Date TT en siecle julien depuis J2000 (reel double).
!
!     t0 : Date du debut de l'intervalle en siecle julien depuis J2000
!          (reel double).
!     dt : Taille de l'intervalle en siecle julien (reel double).
!     nc : Nombre de coefficients (entier, 3<=nc<=10).
!
! --- Common en entree : pol -------------------------------------------
!
!     Tables des coefficients des polynomes Tchebychev (reel double).
!     Le nombre maximum de coefficients par variables est : 10.
!
!     1/ Coefficients des coordonnees rectangulaires equatoriales X,Y,Z
!        de la Lune et du Soleil (reel double) :
!        cal(10,1:3) : X,Y,Z de la Lune en ua.
!        cas(10,1:3) : X,Y,Z du Soleil en ua.
!
!     2/ Coefficients des elements de Bessel.
!        Coordonnees horaires de l'axe passant par le centre de la Terre
!        et parallele a l'axe du cone d'ombre.
!        celb(10,1) : Angle Horaire en radian.
!        celb(10,2) : Declinaison en radian.
!        Coefficients des coordonnees rectangulaires de la Lune dans le
!        systeme de coordonnees fondamental de l'eclipse.
!        celb(10,3) : X en rayon terrestre.
!        celb(10,4) : Y en rayon terrestre.
!        celb(10,5) : Z en rayon terrestre.
!        Coefficients des demi-angles au sommet des cones.
!        celb(10,6) : Demi-angle du cone de penombre en radian.
!        celb(10,7) : Demi-angle du cone d'ombre en radian.
!        Coefficients des rayons des Sections des cones avec le plan
!        fondamental de l'eclipse.
!        celb(10,8) : Rayon du cone de penombre en rayon terrestre.
!        celb(10,9) : Rayon du cone d'ombre en rayon terrestre.
!
! --- Sortie -----------------------------------------------------------
!
!     v1 : Valeurs 1 apres substitution (reel double).
!          Si iv=1, v1(1) : X  Lune en ua.
!                   v1(2) : Y  Lune en ua.
!                   v1(3) : Z  Lune en ua.
!          Si iv=2, v1(1) : H  en radian.
!                   v1(2) : d  en radian.
!                   v1(3) : x  en rayon terrestre.
!                   v1(4) : y  en rayon terrestre.
!                   v1(5) : z  en rayon terrestre.
!                   v1(6) : f1 en radian.
!                   v1(7) : f2 en radian.
!                   v1(8) : u1 en rayon terrestre.
!                   v1(9) : u2 en rayon terrestre.
!
!     v2 : Valeurs 2 apres substitution (reel double).
!          Si iv=1, v2(1) : X  Soleil en ua.
!                   v2(2) : Y  Soleil en ua.
!                   v2(3) : Z  Soleil en ua.
!          Si iv=2, v2(1) : H' en radian / heure.
!                   v2(2) : d' en radian / heure.
!                   v2(3) : x' en rayon terrestre / heure.
!                   v2(4) : y' en rayon terrestre / heure.
!                   v2(5) : z' en rayon terrestre / heure.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension tch(10),tch2(10)
      dimension v1(9),v2(9)
!
      common/pol/
     .cal(10,3),cas(10,3),celb(10,9)
!
! --- Verification du nombre de coefficients ---------------------------
!
      if (nc.lt.3.or.nc.gt.10) stop ' Erreur Polynomes de Tchebychev'
!
      do i=1,9
         v1(i)=0.d0
         v2(i)=0.d0
      enddo
!
! --- Polynomes de Tchebychev ------------------------------------------
!
      x=-1.d0+2.d0*(t-t0)/dt
!
      if (abs(x).gt.1.d0) stop ' Erreur Polynomes de Tchebychev'
!
      tch(1)=1.d0
      tch(2)=x
!
      do i=3,nc
         tch(i)=2.d0*x*tch(i-1)-tch(i-2)
      enddo
!
! --- Selection des coordonnees a calculer -----------------------------
!
      select case (iv)
!
! --- Coordonnees equatoriales de la Lune et du Soleil -----------------
!
      case (1)
!
      do k=1,3
         a1=0.d0
         a2=0.d0
         do i=1,nc
            a1=a1+cal(i,k)*tch(i)
            a2=a2+cas(i,k)*tch(i)
         enddo
         v1(k)=a1
         v2(k)=a2
      enddo
      v1(4)=0.d0
      v2(4)=0.d0
!
! --- Elements de Bessel -----------------------------------------------
!
      case (2)
!
      tch2(1)=0.d0
      tch2(2)=1.d0
!
      do i=3,nc
         tch2(i)=2.d0*x*tch2(i-1)+2.d0*tch(i-1)-tch2(i-2)
      enddo
!
      do j=1,9
         a1=0.d0
         do i=1,nc
            a1=a1+celb(i,j)*tch(i)
         enddo
         v1(j)=a1
      enddo
!
      do j=1,5
         a2=0.d0
         do i=1,nc
            a2=a2+celb(i,j)*tch2(i)
         enddo
         v2(j)=a2*2.d0/dt/36525.d0/24.d0
      enddo
!
      endselect
!
      return
      end
!
!
!
      subroutine LIGNE (te,val,vas,glon,glat,rayc,haut,ierr)
!--------------------------------------------------------------------- -
!
!     Ref : JMG 9904
!
!     Determination d'un point de la ligne de centralite de l'eclipse
!     centrale a partir des coordonnees rectangulaires equatoriales
!     apparentes de la Lune et du Soleil.
!
! --- Entree -----------------------------------------------------------
!
!     te  : Date TT en siecle julien depuis J2000 (reel double).
!
!     val : Coordonnees rectangulaires de la Lune en ua (reel double).
!     vas : Coordonnees rectangulaires du Soleil en ua (reel double).
!
! --- Sortie -----------------------------------------------------------
!
!     glon : Longitude du lieu en degres (reel double).
!     glat : Latitude du lieu en degres (reel double).
!     rayc : Rayon du cone d'ombre sur le lieu en km (reel double).
!     haut : Hauteur du Soleil en degres (reel double).
!
!     ierr : Indice erreur (entier).
!            ierr=0 : Eclipse centrale trouvee.
!            ierr=1 : Eclipse centrale non trouvee.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension val(3),vas(3),xls(3),cp(3)
!
! --- Constantes -------------------------------------------------------
!
      parameter (d2000=2451545.d0,siecle=36525.d0)
      parameter (ua=1.49597870d8)
      parameter (pi=3.141592653589793d0,dpi=2.d0*pi,dgrad=pi/180.d0)
!
! --- Rayons et Applatissement -----------------------------------------
!
      parameter (rl=1738.092d0/ua,rs=696000.d0/ua,rt=6378.140d0/ua)
      parameter (rsl=rs-rl,rt2=rt*rt)
      parameter (uma=1.d0-1.d0/298.257d0,rtp=1.d0/uma,rtp2=rtp*rtp)
!
! --- Recherche du point de centralite ---------------------------------
!
      do i=1,3
         xls(i)=vas(i)-val(i)
      enddo
!
      ca=rtp2*xls(3)*xls(3)
      cb=rtp2*val(3)*xls(3)
      cc=rtp2*val(3)*val(3)
!
      do k=1,2
         ca=ca+xls(k)*xls(k)
         cb=cb+val(k)*xls(k)
         cc=cc+val(k)*val(k)
      enddo
!
      b=cb/ca
      c=(cc-rt2)/ca
      discr=b*b-c
!
      if (discr.lt.0.d0) then
         ierr=1
         return
      endif
!
! --- Coordonnees du point de centralite -------------------------------
!
      beta=-b+sqrt(discr)
      do k=1,3
         cp(k)=xls(k)*beta+val(k)
      enddo
!
! --- Longitude du lieu ------------------------------------------------
!
      alpha=atan2(cp(2),cp(1))
      call TSIDV (te,ts)
      h=ts-alpha
      hp=aint(h/dpi)-1.d0
      gl=h-hp*dpi
      if (gl.gt.pi) gl=gl-dpi
      glon=gl/dgrad
!
! --- Latitude du lieu -------------------------------------------------
!
      tgphi=rtp2*cp(3)/sqrt(cp(1)*cp(1)+cp(2)*cp(2))
      phi=atan(tgphi)
      glat=phi/dgrad
!
! --- Rayon du cone d'ombre sur le lieu --------------------------------
!
      dsl=0.d0
      do k=1,3
         dsl=dsl+xls(k)*xls(k)
      enddo
      sinf=rsl/sqrt(dsl)
      f=asin(sinf)
      r=(rl+beta*rsl)/cos(f)
      rayc=r*ua
!
! --- Hauteur du Soleil (sans refraction) ------------------------------
!
      alpha=atan2(vas(2),vas(1))
      ah=ts-alpha-gl
      dts=sqrt(vas(1)*vas(1)+vas(2)*vas(2)+vas(3)*vas(3))
      delta=asin(vas(3)/dts)
      sinh=sin(phi)*sin(delta)+cos(phi)*cos(delta)*cos(ah)
      haut=asin(sinh)/dgrad
!
      ierr=0
      return
!
      end
!
!
!
      subroutine OMBRE  (v1,v2,rayc,duree,bande)
!-----------------------------------------------------------------------
!
!     Ref : JMG 9904
!
!     Determination de la duree du passage de l'ombre et de la largeur
!     de la bande d'ombre sur un point de la ligne de centralite.
!
!     Les algorithmes utilises sont ceux de J. Meeus.
!     Pour le calcul de la largeur de la bande d'ombre, ils donnent une
!     bonne approximation quand le Soleil est suffisamment haut.
!
! --- Entree -----------------------------------------------------------
!
!     Elements de Bessel :
!     Coordonnees horaires de l'axe passant par le centre de la Terre
!     et parallele a l'axe du cone d'ombre (H,d).
!     Coordonnees de la Lune dans le systeme de coordonnees fondamental
!     de l'eclipse (X,Y,Z).
!     Demi-angles au sommet des cones de penombre et d'ombre (f1,f2).
!     Rayons des Sections des cones de penombre et d'ombre avec le plan
!     fondamental de l'eclipse (u1,u2)
!
!     v1 : Valeurs des elements de Bessel (reel double).
!          v1(1) : H  en radian.
!          v1(2) : d  en radian.
!          v1(3) : x  en rayon terrestre.
!          v1(4) : y  en rayon terrestre.
!          v1(5) : z  en rayon terrestre.
!          v1(6) : f1 en radian.
!          v1(7) : f2 en radian.
!          v1(8) : u1 en rayon terrestre.
!          v1(9) : u2 en rayon terrestre.
!
!     v2 : Valeurs des vitesses des elements de Bessel (reel double).
!          v2(1) : H' en radian / heure.
!          v2(2) : d' en radian / heure.
!          v2(3) : x' en rayon terrestre / heure.
!          v2(4) : y' en rayon terrestre / heure.
!          v2(5) : z' en rayon terrestre / heure.
!
! --- Sortie -----------------------------------------------------------
!
!     duree : Duree du passage de l'ombre en minutes (reel double).
!     bande : Largeur de la bande d'ombre en km (reel double).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      dimension v1(9),v2(9)
!
! --- Initialisation ---------------------------------------------------
!
      duree=1.d9
      bande=1.d9
!
      p=v2(1)
      d=v1(2)
      sd=sin(d)
      cd=cos(d)
!
      x=v1(3)
      y=v1(4)
      xp=v2(3)
      yp=v2(4)
!
      f2=abs(v1(7))
      tanf2=tan(f2)
      u2=abs(v1(9))
      if (rayc.gt.0.d0) u2=-u2
!
      w=1.d0/sqrt(1.d0-0.006694385d0*cd*cd)
      b=yp-p*x*sd
      c=xp+p*y*sd
      y1=w*y
      beta2=1.d0-x*x-y1*y1
      if (beta2.lt.0.d0) return
      beta=sqrt(beta2)
!
! --- Duree du passage de l'ombre --------------------------------------
!
      u2p=u2-beta*tanf2
      a=c-p*beta*cd
      zn2=a*a+b*b
      if (zn2.lt.0.d0) return
      zn=sqrt(zn2)
      duree=120.d0*abs(u2p)/zn
!
! --- Largeur de la bande d'ombre --------------------------------------
!
      z=x*a+y*b
      zk2=beta2+z*z/zn2
      zk=sqrt(zk2)
      bande=12756.d0*abs(u2p)/zk
!
      return
      end
!
!
!
      subroutine EQUAPP (long,lat,r,psi,es,ec,xyz)
! ----------------------------------------------------------------------
!
!     Ref : JMG 9904
!
!     Calcul des coordonnees rectangulaires equatoriales apparentes de
!     la Lune et du Soleil a partir des coordonnees ecliptiques moyennes
!     corrigees de l'aberration (longitude, latitude et distance).
!
! --- Entree -----------------------------------------------------------
!
!     long  : Longitude en radian (reel double).
!     lat   : Latitude en radian (reel double).
!     r     : Distance a la Terre en ua (reel double).
!
!     psi   : Nutation en longitude en radian (reel double).
!     es    : Sinus de l'obliquite vraie (reel double).
!     ec    : Cosinus de l'obliquite vraie (reel double).
!
! --- Sortie -----------------------------------------------------------
!
!     xyz   : Coordonnees rectangulaires equatoriales apparentes.
!             xyz(1) : X en ua (reel double).
!             xyz(2) : Y en ua (reel double).
!             xyz(3) : Z en ua (reel double).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      double precision long,lat
!
      dimension xyz(3)
!
! --- Longitude et Latitude apparentes ---------------------------------
!
      a=long+psi
      ca=cos(a)
      sa=sin(a)
!
      if ((abs(lat)-1.d-10).lt.0.d0) then
         sb=0.d0
         cb=r
      else
         sb=r*sin(lat)
         cb=r*cos(lat)
      endif
!
! --- Coordonnees rectangulaires equatoriales apparentes ---------------
!
      p=cb*sa
      x=cb*ca
      y=p*ec-sb*es
      z=p*es+sb*ec
!
      xyz(1)=x
      xyz(2)=y
      xyz(3)=z
!
      return
      end
!
!
!
      subroutine BESSEL (val,vas,tsv,velb)
!--------------------------------------------------------------------- -
!
!     Ref : JMG 9904
!
!     Determination des elements de Bessel.
!
! --- Entree -----------------------------------------------------------
!
!     val  : Coordonnees rectangulaires equatoriales apparentes de la
!            Lune (reel double).
!            val(1) : X en ua.
!            val(2) : Y an ua.
!            val(3) : Z en ua.
!
!     vas  : Coordonnees rectangulaires equatoriales apparentes du
!            Soleil (reel double).
!            vas(1) : X en ua.
!            vas(2) : Y an ua.
!            vas(3) : Z en ua.
!
!     tsv   : Temps Sideral Vrai en radian (reel double).
!
! --- Sortie -----------------------------------------------------------
!
!     velb  : Elements de Bessel (reel double).
!
!             Coordonnees horaires de l'axe passant par le centre de la
!             de la Terre et parallele a l'axe du cone d'ombre.
!             velb(1) : H  Angle Horaire en radian.
!             velb(2) : d  Declinaison en radian.
!             Coordonnees rectangulaires de la Lune dans systeme de
!             coordonnees fondamental de l'eclipse.
!             velb(3) : X  en rayon terrestre.
!             velb(4) : Y  en rayon terrestre.
!             velb(5) : Z  en rayon terrestre.
!             Demi-angles au sommet des cones de penombre et d'ombre.
!             velb(6) : f1 Demi-angle du cone de penombre en radian.
!             velb(7) : f2 Demi-angle du cone d'ombre en radian.
!             Rayons des Sections des cones de penombre et d'ombre avec
!             le plan fondamental de l'eclipse.
!             velb(8) : u1 Rayon du cone de penombre en rayon terrestre.
!             velb(9) : u2 Rayon du cone d'ombre en rayon terrestre.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension val(3),vas(3),velb(9)
!
! --- Constantes -------------------------------------------------------
!
      parameter (ua=1.49597870d8)
      parameter (pi=3.141592653589793d0,dpi=2.d0*pi)
      parameter (sdrad=pi/180.d0/3600.d0)
      parameter (pi0=8.794148d0*sdrad,s0=(15.d0*60.d0+59.63d0)*sdrad)
      parameter (rl=1738.092d0,rt=6378.140d0,rlt=rl/rt)
!
! --- Initialisation ---------------------------------------------------
!
      dsl=0.d0
      dtl=0.d0
      dts=0.d0
!
      do i=1,3
         u=vas(i)-val(i)
         dsl=dsl+u*u
         dtl=dtl+val(i)*val(i)
         dts=dts+vas(i)*vas(i)
      enddo
!
      dsl=sqrt(dsl)
      dtl=sqrt(dtl)
      dts=sqrt(dts)
!
      al=atan2(val(2),val(1))
      dl=asin(val(3)/dtl)
      as=atan2(vas(2),vas(1))
      ds=asin(vas(3)/dts)
!
      g=dsl/dts
      h=dtl/dts
!
      cal=cos(al)
      sal=sin(al)
      cas=cos(as)
      sas=sin(as)
      cdl=cos(dl)
      sdl=sin(dl)
      cds=cos(ds)
      sds=sin(ds)
!
! --- Coordonnees de l'axe du cone d'ombre -----------------------------
!
      asc=atan2((cds*sas-h*cdl*sal),(cds*cas-h*cdl*cal))
      ah=tsv-asc
      ah=mod(ah,dpi)
      dec=asin((sds-h*sdl)/g)
!
      velb(1)=ah
      velb(2)=dec
!
! --- Coordonnees rectangulaires de la Lune ----------------------------
!
      cd=cos(dec)
      sd=sin(dec)
      ca=cos(al-asc)
      sa=sin(al-asc)
!
      x=dtl*cdl*sa*ua/rt
      y=dtl*(sdl*cd-cdl*sd*ca)*ua/rt
      z=dtl*(sdl*sd+cdl*cd*ca)*ua/rt
!
      velb(3)=x
      velb(4)=y
      velb(5)=z
!
! --- Angles et rayons des cones ---------------------------------------
!
      sinf1=(sin(s0)+rlt*sin(pi0))/g/dts
      sinf2=(sin(s0)-rlt*sin(pi0))/g/dts
!
      f1=asin(sinf1)
      f2=asin(sinf2)
      if (f1.lt.0.d0) f1=-f1
      if (f2.gt.0.d0) f2=-f2
!
      u1=z*tan(f1)+rlt/cos(f1)
      u2=z*tan(f2)+rlt/cos(f2)
!
      velb(6)=f1
      velb(7)=f2
      velb(8)=u1
      velb(9)=u2
!
      return
      end
!
!
!
      subroutine NUTAT (t,psi,es,ec)
! ----------------------------------------------------------------------
!
!     Ref : JMG 9904
!
!     Calcul de la nutation en longitude et obliquite.
!
!     On utilise les formules simplifiees de P. Bretagnon et J.L. Simon.
!
! --- Entree -----------------------------------------------------------
!
!     t   : Date TT en siecle julien depuis J2000 (reel double).
!
! --- Sortie -----------------------------------------------------------
!
!     psi : Nutation en longitude en radian (reel double).
!     es  : Sinus de l'obliquite vraie (reel double).
!     ec  : Cosinus de l'obliquite vraie (reel double).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension ie(5)
!
! --- Initialisation ---------------------------------------------------
!
      data ie/-226938,-75,96926,-2491,-12104/
!
      u=t/100.d0
!
      a1=2.18d0+(0.36d0*u-3375.70d0)*u
      a2=3.51d0+(0.10*u+125666.39d0)*u
      s1=sin(a1)
      s2=sin(a2)
      c1=cos(a1)
      c2=cos(a2)
!
! --- Nutation en longitude --------------------------------------------
!
      psi=-834.d-7*s1-64.d-7*s2
!
      el=0.d0
      do i=5,1,-1
         el=u*(el+ie(i))
      enddo
!
! --- Obliquite vraie --------------------------------------------------
!
      epsi=0.4090928d0+1.d-7*(el+446*c1+28*c2)
!
      es=sin(epsi)
      ec=cos(epsi)
!
      return
      end
!
!
!
      subroutine TSIDV (t,ts)
! ----------------------------------------------------------------------
!
!     Ref : JMG 9904
!
!     Calcul du Temps Sideral vrai au meridien origine.
!
! --- Entree -----------------------------------------------------------
!
!     t  : Date en siecle julien depuis J2000 (reel double).
!
! --- Sortie -----------------------------------------------------------
!
!     ts : Temps Sideral vrai en radian (reel double).
!
! --- Remarque ---------------------------------------------------------
!
!     Cette subroutine donne le Temps Sideral vrai a Greenwich si la
!     date t est exprimee dans l'echelle du Temps Universel (UT).
!     Dans le cas ou le temps t est exprime en Temps Terrestre (TT), il
!     s'agit du Temps sideral vrai pour le "meridien des Ephemerides".
!     La longitude L du meridien de Greenwich par rapport a celui des
!     Ephemerides est : L=(TT-UT)*1.0027379, a l'Ouest si L est positif
!     et a l'Est si L est negatif.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
! --- Parametres -------------------------------------------------------
!
      parameter (siecle=3600.d0*24.d0*36525.d0)
      parameter (tstm=1.0027379d0*siecle)
      parameter (shrad=15.d0/206264.8062470964d0)
!
! --- Calcul du Temps Sideral moyen au meridien origine ----------------
!
      dj=aint(t*36525.d0)
      t0=(dj+0.5d0)/36525.d0
      ts0=24110.54841d0+t0*(8640184.812866d0+(0.093104d0-6.2d-6*t0)*t0)
      tsmg=ts0+(t-t0)*tstm
      tsmg=tsmg*shrad
!
! --- Calcul du Temps Sideral vrai -------------------------------------
!
      call NUTAT (t,psi,es,ec)
      ts=tsmg+psi*ec
!
      return
      end
!
!
!
      subroutine DTTUT (tdj,dift)
!-----------------------------------------------------------------------
!
!     Ref : JMG 9904
!
!     Difference entre Temps Terrestre (TT) et Temps Universel (UT).
!
! --- Entree -----------------------------------------------------------
!
!     tdj  : Date julienne (reel dp).
!
! --- Sortie -----------------------------------------------------------
!
!     dift : Valeur de la difference TT-UT en secondes (reel dp).
!
! --- Remarque ---------------------------------------------------------
!
!     Pour les annees anterieures a 1955, on utilise les formules et
!     tables de Stephenson et Morrison (1984) corrigees de la valeur
!     de l'acceleration de la Lune due aux marees -25.738"/siecle.
!     Entre 1955 et 2000, on donne des valeurs deduites a la seconde
!     pres de celles fournies par l'IERS pour TAI-UTC.
!     Au-dela de 2000, on fait une estimation empirique calculee a
!     partir des valeurs observees sur l'intervalle 1800-2000.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension ns1(37),ns2(175),ns3(47)
!
      parameter (d2000=2451545.d0,siecle=36525.d0)
      parameter (t948=-10.519658d0)
      parameter (t1600=-3.999932d0)
      parameter (t1780=-2.199996d0)
      parameter (t1955=-0.450007d0)
      parameter (t2001=+0.010021d0)
!
! --- Initialisation ---------------------------------------------------
!
      data ns1/
     .102,99,96,94,91,88,85,72,62,54,48,43,37,32,26,21,16,12,10,09,
     . 09,09,10,10,11,11,11,12,12,13,13,14,15,16,16,17,17/
!
      data ns2/
     .17,17,17,17,17,17,17,17,17,17,17,17,17,16,16,16,15,15,14,14,
     .14,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,12,12,12,
     .12,12,11,11,11,10,10,09,09,08,08,07,07,06,06,06,06,06,06,06,
     .06,06,06,06,06,06,07,07,07,07,07,07,07,07,08,08,08,08,08,08,
     .08,08,08,07,06,06,05,04,03,02,02,00,-1,-1,-3,-3,-4,-5,-5,-5,
     .-5,-5,-5,-5,-5,-6,-6,-6,-6,-6,-6,-6,-6,-7,-6,-6,-6,-6,-5,-4,
     .-3,-2,00,01,03,04,05,06,08,09,10,12,13,15,16,17,18,19,20,21,
     .21,22,22,23,23,24,24,24,24,24,24,24,24,24,24,24,24,24,24,24,
     .24,25,25,26,26,27,27,28,28,29,29,30,30,30,31/
!
      data ns3/
     .31,31,32,32,33,33,34,34,34,35,36,37,37,38,39,40,41,42,43,44,
     .45,46,48,49,50,51,51,52,53,54,54,55,55,56,56,57,58,58,59,60,
     .61,62,62,63,64,65,66/
!
! --- Calcul -----------------------------------------------------------
!
      t=(tdj-d2000)/siecle
!
      if (t.lt.t948) then
         dift=2177.d0+t*(497.d0+t*44.1d0)
         return
      endif
!
      if (t.lt.t1600) then
         dift=102.d0+t*(102.d0+t*25.3d0)
         return
      endif
!
      if (t.lt.t1780) then
         i=aint((t-t1600)/0.05d0)+1
         t0=t1600+(i-1)*0.05d0
         x=ns1(i)+(t-t0)*(ns1(i+1)-ns1(i))/0.05d0
         dift=x-0.2d0*t*t
         return
      endif
!
      if (t.lt.t1955) then
         i=aint((t-t1780)/0.01d0)+1
         dift=ns2(i)-0.2d0*t*t
         return
      endif
!
      if (t.lt.t2001) then
         i=aint((t-t1955)/0.01d0)+1
         dift=ns3(i)
         return
      endif
!
      dift=65.d0+t*(91.d0+t*32.8d0)
!
      return
      end
!
!
!
      subroutine DJDATE (tdj,icode,jour,mois,kan,nh,nm,ns)
!-----------------------------------------------------------------------
!
!     Ref : JMG 9904
!
!     Conversion Date julienne -> Date du calendrier (Gregorien/Julien).
!
!     On utilise le formulaire de J. Meeus.
!
! --- Entree ------------------------------------------------------------
!
!     tdj    :  Date julienne (reel dp).
!
!     icode  :  Indice code de la conversion (entier)
!               icode=0 conversion tdj -> jour-mois-kan et 00-00-00.
!               icode=1 conversion tdj -> jour-mois-kan et nh-00-00.
!               icode=2 conversion tdj -> jour-mois-kan et nh-nm-00.
!               icode=3 conversion tdj -> jour-mois-kan et nh-nm-ns.
!
! --- Sortie -----------------------------------------------------------
!
!     jour   :  Quantieme du mois (entier).
!     mois   :  Numero du mois (entier).
!     kan    :  Millesime astronomique de l'annee (entier).
!
!     nh     :  Nombre d'heures (entier).
!     nm     :  Nombre de minutes (entier).
!     ns     :  Nombre de secondes (entier).
!
! --- Remarque ---------------------------------------------------------
!
!     Date du calendier (Gregorien/Julien) :
!     Le calendier Gregorien commence le 15/10/1582 a 0h (DJ2299160.5).
!     Avant cette date on utilise le calendrier Julien.
!
!     La date julienne tdj doit etre superieure ou egale a 0.
!     Si tdj<0, alors jour=mois=kan=nh=nm=ns=0.
!
!     L'annee kan est comptee en millesimes astronomiques.
!     L'annee 1 avant J.C. des historiens correspond a l'an 0 des
!     astronomes et ainsi de suite.
!     Ex : l'an 80 avant J.C. correspond a l'an astronomique -79.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
! --- Initialisation ---------------------------------------------------
!
      jour=0
      mois=0
      kan=0
      nh=0
      nm=0
      ns=0
!
      if (tdj.lt.0.d0.or.icode.lt.0.or.icode.gt.3) return
      if (icode.eq.0) dt=0.0d0
      if (icode.eq.1) dt=0.5d0/24.d0
      if (icode.eq.2) dt=0.5d0/1440.d0
      if (icode.eq.3) dt=0.5d0/86400.d0
!
! --- Conversion -------------------------------------------------------
!
      t=tdj+dt+0.5d0
      z=int(t)
      f=t-z
!
      if (z.lt.2299161.d0) then
         a=z
      else
         alpha=int((z-1867216.25d0)/36524.25d0)
         a=z+1.d0+alpha-int(alpha/4.d0)
      endif
!
      b=a+1524.d0
      c=int((b-122.1d0)/365.25d0)
      d=int(365.25d0*c)
      e=int((b-d)/30.6001d0)
!
      jour=b-d-int(30.6001d0*e)
!
      if (e.lt.13.5d0) then
         mois=e-1.d0
      else
         mois=e-13.d0
      endif
!
      if (mois.lt.2.5d0) then
         kan=c-4715.d0
      else
         kan=c-4716.d0
      endif
!
      f=f*24.d0
      nh=f
      f=(f-nh)*60.d0
      nm=f
      f=(f-nm)*60.d0
      ns=f
!
      return
      end
!
!
!
      subroutine DATEDJ (jour,mois,kan,nh,nm,ns,tdj)
!-----------------------------------------------------------------------
!
!     Ref : JMG 9904
!
!     Conversion Date du calendrier (Gregorien/Julien) -> Date julienne.
!
!     On utilise le formulaire de J. Meeus.
!
! --- Entree ------------------------------------------------------------
!
!     jour   :  Quantieme du mois (entier).
!     mois   :  Numero du mois (entier).
!     kan    :  Millesime astronomique de l'annee (entier).
!
!     nh     :  Nombre d'heures (entier).
!     nm     :  Nombre de minutes (entier).
!     ns     :  Nombre de secondes (entier).
!
! --- Sortie -----------------------------------------------------------
!
!     tdj     : Date julienne (reel dp).
!
! --- Remarque ---------------------------------------------------------
!
!     Date du calendier (Gregorien/Julien) :
!     Le calendier Gregorien commence le 15/10/1582 a 0h (DJ2299160.5).
!     Avant cette date on utilise le calendrier Julien.
!
!     Cette fonction n'est valable que pour les dates posterieures au
!     1/1/4713 12h avant J.C.
!
!     L'annee kan est comptee en millesimes astronomiques.
!     L'annee 1 avant J.C. des historiens correspond a l'an 0 des
!     astronomes et ainsi de suite.
!     Ex : l'an 80 avant J.C. correspond a l'an astronomique -79.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
! --- Initialisation ---------------------------------------------------
!
      tdj=0.d0
!
      if ((kan.lt.-4713)
     .   .or.(mois.lt.1.or.mois.gt.12)
     .   .or.(jour.lt.1.or.jour.gt.31)
     .   .or.(nh.lt.0.or.nh.gt.24)
     .   .or.(nm.lt.0.or.nm.gt.59)
     .   .or.(ns.lt.0.or.ns.gt.59)) return
!
      a=0.d0
      b=0.d0
      c=0.d0
!
! --- Conversion -------------------------------------------------------
!
      if (mois.gt.2) then
         y=kan
         m=mois
      else
         y=kan-1.d0
         m=mois+12
      endif
!
      if (y.lt.0) then
         c=-0.75d0
      else
         idate=kan*10000+mois*100+jour
         if (idate.ge.15821015) then
            a=int(y/100.d0)
            b=2.d0-a+int(a/4.d0)
         endif
      endif
!
      tdj=int(365.25d0*y+c)
     .   +int(30.6001d0*(m+1))
     .   +jour+b
     .   +dfloat(nh)/24.d0
     .   +dfloat(nm)/1440.d0
     .   +dfloat(ns)/86400.d0
     .   +1720994.5d0
!
      return
      end
!
!
!
      subroutine CONTV (t,n,v)
!-----------------------------------------------------------------------
!
!     Ref : GF9803
!
!     Rendre continue une suite de valeurs exprimees en radians autour
!     d'une valeur moyenne v.
!
! --- Entree -----------------------------------------------------------
!
!     t  : Table des valeurs en radians (reel dp)
!     n  : Nombre de valeurs (entier).
!     v  : Valeur moyenne en radians (reel dp).
!
! --- Sortie -----------------------------------------------------------
!
!     t  : Table des valeurs en radians (reel dp)
!
!---- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension t(n)
!
      parameter (pi=3.141592653589793d0,dpi=2.d0*pi)
!
!---- Calcul -----------------------------------------------------------
!
      d=t(1)-v
      t(1)=t(1)-dint(d/dpi)*dpi
      d=t(1)-v
      if (abs(d).gt.pi) t(1)=t(1)-dsign(1.d0,d)*dpi
!
      do i=2,n
         d=t(i)-t(i-1)
         do while (abs(d).gt.pi)
            t(i)=t(i)-dpi*d/abs(d)
            d=t(i)-t(i-1)
         enddo
      enddo
!
      return
      end
!
!
!
      subroutine EFFECR
!-----------------------------------------------------------------------
!
!     Ref : JMG 9904
!
!     Effacement de l'ecran.
!
! --- Effacement de l'ecran sous DOS -----------------------------------
!
      call system ('cls')
!
! --- Effacement de l'ecran sous UNIX ----------------------------------
!**
!**   call system ('clear')
!**
!
      return
      end
